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Disclaimers 


The  findings  In  this  report  are  not  to  be  construed  as  an 
official  Department  of  the  Army  position,  unless  so  desig¬ 
nated  by  other  authorized  documents. 

The  citation  of  trade  names  and  names  of  manufacturers  in 
this  report  Is  not  to  be  construed  as  official  Government 
indorsement  or  approval  of  commercial  products  or  services 
referenced  herein. 


Disposition 

Destroy  this  report  when  it  is  no  longer  needed.  Do  not 
return  it  to  the  originator. 
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1.  INTRODUCTION 


A  continued  fraction  technique1  2  has  been  developed  for  the  calculation  of 
the  ratios  of  spherical  Bessel  functions  of  complex  argument  used  In  Mle 
calculations.3  4  When  the  spherical  Bessel  functions  of  complex  argument 
themselves  are  required,  the  same  technique  may  be  applied  with  several 
simplifications.  Spherical  Bessel  functions  of  complex  argument  have  wide 
ranging  applications  in  physics,  electronics,  and  meteorology.  The  purpose  of 
this  report  is  to  develop  the  simplified  formulae  for  the  calculation  of 
spherical  Bessel  functions  of  complex  argument  and  to  present  several 
algorithms  that  correctly  calculate  the  functions.  The  technique  is 
essentially  based  on  a  new  method  of  evaluating  continued  fractions. 

There  are  now  four  distinctly  different  methods  that  allow  one  to  evaluate  a 
continued  fraction.  Three  of  these  methods  are  detailed  in  Gautschi,5  and  the 
fourth  in  Lentz.1  The  Lentz6  continued  fraction  algorithm  Is  of  much  more 
general  application  than  just  spherical  Bessel  functions.  This  algorithm 
essentially  transforms  a  continued  fraction  into  an  infinite  product,  whereas 
the  previous  three  methods  convert  to  infinite  series  or  recursions  that  have 
stability  problems. 


2.  SIMPLE  CONTINUED  FRACTION 

A  simple  continued  fraction  may  be  defined  as: 


1 W.  J.  Lentz,  1973,  A  New  Method  of  Confuting  Spherical  Bessel  Functions  of 
Complex  Argument  with  Tables,  ECUR  5509,  AD  /b/2Z3,  US  Army  Atmospheric 
Sciences  Laboratory,  White  Sands  Missile  Range,  NM 

2 W .  J.  Lentz,  1976,  "Generating  Bessel  Functions  In  Mie  Scattering  Calcula¬ 
tions  Using  Continued  Fractions,"  Appl  Opt,  15:668 

3W.  J.  Wlscombe,  1980,  'Unproved  Mle  Scattering  Algorithms,"  Appl  Opt,  19:1505 

HG.  Grehan  and  G.  Gouesbet,  1979,  "Mle  Theory  Calculations:  New  Progress, 
with  Emphasis  on  Particle  Sizing,"  Appl  Opt,  18:3489 

'JW.  Gautschi,  1967,  SIAM,  Rev  9,  No  1 


bW.  J.  Lentz,  1982,  A  Si mpl i f 1 cati on  of  Lentz ' s  A1 gori thm ,  ASL-TR-0118,  US 
Army  Atmospheric  Sciences  Laboratory.  White  Sands  Missile  Range.  NM  i 
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where  none  of  the  a^  may  be  zero.  The  a^  may  be  complex  or  negative,  and  the 
nth  convergent  of  the  continued  fraction  is  defined  as: 


Inis  simple  continued  fraction  might  better  be  called  a  positive  simple 
continued  fraction  because  all  the  terms  are  added.  One  method  of  eva^ati 

equation  (2)  is  to  start  at  an,  take  the  reciprocal  and  add  to  an  ^  until  a( 
is  reached.  The  only  trick  is  to  know  at  which  an  to  start  for  a  given 

accuracy.  Unfortunately,  two  terms  in  the  iteration  process  may  cancel,  or 
nearly  cancel,  causing  unpredictable  errors  in  the  final  answer.  Let  us 
therefore  review  a  new  technique  to  evaluate  continued  fractions.2 

set  us  make  a  further  change  of  notation  to  define  a  simple  continued  fraction 
■no re  compactly. 

Fn  =  [al»  a2.  a3 . an]  (3J 

Lentz2  has  shown  that  this  continued  fraction  may  be  transformed  to  an 
equivalent  form,  which  begins  at  a,  rather  than  some  indeterminate  an. 


Although  this  formula  may  appear  complicated,  one  generates  successive 
numerators  and  denominators  by  taking  the  reciprocal  and  adding  to  the 
previous  numerator  or  denominator.  A  numerical  example  and  an  algorithm  for 
equation  (4)  may  be  found  in  reference  6. 


W.  J.  Lentz,  1976,  "Generating  Bessel  Functions  in  Mie  Scattering 
Calculations  Using  Continued  Fractions,"  Appl  Opt,  15:668 

’  w.  J.  Lentz,  1982,  A  Simplification  of  Lentz's  Algorithm,  ASL-TR-0118,  US 
Army  Atmospheric  Sciences  Laboratory,  White  Sands  Missile  ftange,  NM 
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Since  the  a^  might  be  negative,  one  will  note  that  any  given  term 

[a..,  a^  ^ . al(or  2)^  be  inaccurate  due  to  cancellation  of 

common  digits  in  the  sum.  In  this  case,  an  algorithm  improvement  has  been 

devi sed.2 
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'  ’S  tne  product  of  the  next  two  terms.  If  B  is  identically  zero,  then  ;  is 
one  2  the  next  term  is  then: 
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Note  that  the  above  formulae  apply  to  both  the  numerator  and  the  denominator 
since  the  form  is  independent  of  whether  a  or  a^  is  used.  In  fact,  accuracy 
can  be  increased  by  using  the  algorithm  improvement  any  time  a  common  number 

of  digits  cancel  to  reducing  accuracy  in  the  generation  of  e.  When  a  given 
numerator  such  as  B  is  equal  to  its  correspondi ng  denominator  to  the  number  of 
desired  digits,  the  fraction  may  be  said  to  have  converged  to  the  desired 
accuracy.  Reasonable  accuracy  may  also  be  obtained  by  only  considering  the 
case  when  B  is  identically  zero,  and  the  simplier  algorithm  will  improve 
execution  time  significantly. 


2 W .  J.  Lentz,  1976,  "Generating  Bessel  Functions  in  Mie  Scattering  Calcula¬ 
tions  Using  Continued  Fractions,"  Ap pi  Opt,  15:668 


3.  SIMPLIFIED  FORMULAE 


The  continued  fraction  representation  for  the  ratio  of  Bessel  functions  ot 
consecutive  order  is  given  by: 


(7) 


ere  the  an  =  (-l)n+*2(v+n-l)/x,  v  is  the  order,  and  x  is  the  argument.  A 

simplification  is  possible  in  this  case  because  every  other  term  is  of 
opposite  sign.  The  alternating  sign  can  be  eliminated  by  transforming  the 
simple  continued  fraction  to  a  simple  negative  continued  fraction  using 
equation  (5)  of  3.10.1  of  Abramowitz.7 
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which  will  be  defined  as  a  negative  simple  continued  fraction.  None  of  the  bi 
may  be  zero,  but  they  may  be  complex  or  negative.  In  like  manner,  a  new 
notation  may  be  defined  for  the  nth  convergent  Gn. 


Gn  =  [t>1’  bz’  •  bn3 

=  [blt  -b,.  b3,  -b4 ,  b5 ,  -b6,  ...  ,  bn] 


(10) 


7M.  Abramowitz,  F.  W.  J.  Olver,  and  H.  A.  Antosiewicz,  1965,  Handbook  of 
Mathematical  Functions,  M.  Abramowitz  and  I.  A.  Stegun,  editor?;  Applied 
Mathematics  Series  55,  chapters  3,  9,  and  10,  US  Government  Printing  Office, 
Washington,  OC 
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where  the  first  equation  is  a  negative  simple  continued  fraction  and  the 
second  equation  an  equivalent  positive  simple  continued  fraction. 

The  need  for  this  additional  complexity  will  become  apparent  when  the  formulae 
for  the  spherical  Bessel  functions  is  derived. 

The  equivalent  algorithm  for  negative  continued  fractions  is 
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The  algorithm  improvement  for  negative  continued  fractions  follows  directly 
from  the  definition  of  the  Dm. 


C  s  By  =  bm+1e  -  1  *12) 

l  *  bm+2  -  B/c  (13) 


The  above  formulae  may  apply  to  either  the  numerator  or  denominator  or  both 
completely  independently . 


4.  SPHERICAL  BESSEL  FUNCTIONS 

One  straightforward  way  to  calculate  the  consecutive  orders  of  Bessel 
functions  by  the  continued  fraction  technique  is  to  multiply  appropriate 
ratios  by  an  initial  value 
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The  initial  value  for  spherical  Bessel  functions  of  complex  argument  is  easily 
calculated  from  equation  10.1.11  of  Abramowitz.7 


jQ(x)  =  sin(x)/x  (15) 

In  terms  of  a  negative  continued  fraction  bm  =  2  ( v  +-m- 1 ) /x  where  v  -  n  <-  0.5  i  ■> 
the  order  and  m  is  the  number  of  the  term  in  the  continued  fraction. 


Let  bn  m  =  2(n  +  m  -  0.5)/x 


This  definition  is  exactly  equivalent  to  the  previous  b^,  and  it  allows  one  to 
keep  track  of  both  the  order  of  the  Bessel  function  ratio  and  the  order  of  the 
continued  fraction  term.  Using  this  notation,  the  product  of  n  continued 
fraction  ratios  and  an  initial  value  becomes: 


1/j  -  l/jn  [-]  [b  ][b  b.  j  ...  [o  l[b  ,b  ]  ... 

n  0  1,1  1,2  i,l  n.l  n,2  n,l 


[b.  J  .. 


[b  ][b  b  ] 
n,2  n,3  n,2 


Some  of  the  terms  in  the  numerator  and  denominator  are  identical  and  cancel; 
for  example,  [b1  2]  =  [b2  j]  .  The  terms  that  are  generally  identical  are: 


[bs,m+l3r 


^■bs-l,m^ 


[bs,w2’bs,m+l]niim  =  Cb 


s-1 ,m+l ,bs-l,m^ 


A  simple  rule  to  use  to  understand  the  terms  that  are  identical  is  to  note 
that  lower  order  denominators  cancel  the  next  higher  order  numerators.  All  of 
the  denominators  corresponding  to  the  ratio  of  Bessel  functions  of  a  given 
order  cancel  all  of  the  numerators  of  the  next  higher  order  ratio  of  Bessel 
functions.  All  that  is  left  is  the  numerator  of  the  ratio  J 0/ J  i  and  the 
denominator  of  the  ratio  l /jn .  This  cancellation  is  not  approximate. 
Every  term  is  exactly  identical  If  the  proper  terms  are  considered.  Let  us 
consider  a  numerical  example  for  clarification.  Let  x  =  1.0  in  jg(x). 


7M.  Abramowitz,  F.  W.  J.  Olver,  and  H.  A.  Antosiewicz,  1965,  Handbook  of 
Mathematical  Functions,  M.  Abramowtiz  and  I.  A.  Stegun,  editors,  Applied 
Mathematics  Series  55,  chapters  3,  9,  and  10,  US  Government  Printing  Office, 
Washington,  DC 
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1/jy  =  (3) (4. 666) (6. 785714286) (8. 852631589) (10. 88703924) (12. 90814766) 


( 14.92252955) (16.93298723) (18.94094367) (20. 94720432) (22. 95226093) 
(24.95643131) (26 .95993017)  *  l/sin(l)  =  (1.4913765  x  10"9)'1 


Note  that  there  are  nine  terms  before  the  denominator  begins,  and  convergence 
is  rapid  once  the  continued  fraction  is  begun. 


It  is  important  to  realize  that  the  error  correction  algorithm  improveme  it 
applies  to  the  first  nine  terms  above  as  well  as  to  the  terms  in  the  numerator 
am  denominator  of  the  complete  continued  fraction.  The  application  of  the 
algorithm  improvement  is  completely  independent  of  what  is  happening  else¬ 
where,  and  in  a  properly  constructed  algorithm,  no  loss  of  accuracy  will  occur 
due  to  cancellation  of  digits  in  addition  or  subtraction. 


6.  f-URTRAN  PROGRAMS 

first  FORTRAN  IV  program  m  the  appendix  is  a  simplified  version  of  the 
spherical  Bessel  function  formulae  developed  above.  The  only  error  check  is 
for  a  term  being  identically  zero  in  either  the  numerator  or  denominator.  In 
that  case,  the  error  bypass  is  implemented  through  the  algorithm  improvement 
as  detailed  in  reference  6.  If  a  given  term  is  zero,  the  product  of  that  term 
and  the  next  in  the  algorithm  improvement  is  one.  Therefore,  it  is  sufficient 
to  simply  skip  the  product  and  reciprocal. 

"he  second  FORTRAN  IV  program  is  essentially  the  same  program  modified 
slightly  for  complex  arguments.  The  test  for  convergence  is  also  appro¬ 
priately  modified.  Both  of  these  programs  have  been  checked  against  tables1 
of  spherical  3essel  functions  and  found  to  be  more  accurate  than  the  tables. 
Remember  that  the  final  answer  is  scaled  by  powers  of  10**50  so  that  overflow 
will  not  occur  in  tne  computer.  The  displayed  answer  mist  be  adjusted  to 
account  for  this  if  the  values  are  used  in  confutations. 


bW.  J.  Lentz,  1982,  A  Simplification  of  Lentz's  Algorithm,  ASL-TR-0118,  US 
Army  Atmospheric  Sciences  "laboratory,  dhi  te~ bands  Missile"  Range,  NM 


W.  .).  Lentz,  19/3,  A  New  Method  of  Computing  Spherical  Bessel  Functions  of 
Complex  Argument  wi  tR  Tab  1  es~  E!C0M  5509,  AD  /b it'll ,  OS  Army  Atmospheric 
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APPENDIX 


FORTRAN  ALGOR  IT}**  LISTINGS 


r 


NOVA  FORTRAN  5,  Version  6.16  --  Monday.  Apr  i  1  23,  1994  11:27:53  AM 

FBESSEL.FR 

:  COMPILER  DOUBLE  PRECISION 

:  C  PROGRAM  FBESSEl  CALCULATES  SPHERiCAL  I3FSSEL  FUNCTIONS  OF  THE  FIRST 

:  C  KIND  BY  LENTZ  CONTINUED  FRACTION  UITH  A  CHECK  FOR  ZERO  ONLY 

:  C  IN  THE  NUMERATOR  AND  DENOMINATOR  CALCULATIONS 

:  C  SEE  APPLIED  OPTICS  VOL  15,  NO  3  MARCH  1976  P668. 

INTEGER  ICOUNT.N 

REAL  V, X INC , AN, JO , JN, X.NUM, DEN, PD T, SCALE 

:  1  TYPE  'INPUT  ORDER,  ARGUMENT  OF  SPHERICAL  BESSEL  FUNCTION* 

:  ACCEPT  N.X 

PDT-i.0 

:  V-FLOATCN1+0.5 

NUM-0 . 0 
DEN -0.0 
SCALE *8.0 
:  AN-1.0/X 

;  XlNC«2.e/X 

:  ICOUNT-0 

:  C  IF  N  IS  ZERO.  THE  INITIAL  VALUE  SINCX)XX  IS  USED 

IFCN.EQ.0)  GOTO  20 
:  17  AN-AN+XINC 

NUM-AN-NUM 

IF (NUM. EQ . 0 . 0)  GOTO  13 
:  PDT-PDTwNUM 

:  NUM-1.0/NUM 

IF (ABSIPDT) ,LT. 1 . 0D50)  GOTO  18 
:  PDT-PDT*! .0D-50 

:  SCALE-SCALE-50. 

:  15  ICOUNT- ICOUNT+1 

:  IF ( I COUNT. LE . N)  GOTO  19 

:  DEN-AN-DEN 

:  IF (DEN . EO . 0 . 0)  GOTO  19 

1;  PDT-PDTxDEN 

i:  DEN-1.0XDEN 

19  IF  (NUM.NE .  DEN .  OR,  I  COUNT.  LT.  5)  GOTO  17 

20  JO«SIN(X)/X 

1:  JN-JO/PDT 

I:  URITEC10, 100)  JN. SCALE. X. V. ICOUNT 

I:  100  FORMAT ( 1H  . " JN- " . D20 . 10 . *  POUER  OF  10  SCALE-*. F6. 0.x, 1H  . *X«*. (E14.6) 

:  X  *  V-".F6. 1. "  COUNT-'. 15) 

! :  GOTO  1 
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COMPILER  DOUBLE  PRECISION 

C  PROGRAM  FCBESSEL  CALCULATES  BESSEL  FUNd.ONS  BY  FIRST  RATIO  METHOD 

C  AND  SIMPLIFICATIONS  UITH  ONLY  A  CHECK  FOP  ZERO.  NO  ACCURACY  IMPROVE MEN 

C  IS  IMPLEMENTED 

r  ¥'#'**#*  +  *  +  ■*'**'**'>*;*  **'*#***■  t 

REAL  COUNT. V. SCALE 
COMPLEX  X INC. AN. JO. JN.X.NUM. DEh.PDT 

C  *-■**+  mwX-*W*~*W*+***mXXXor+Xj¥.#*#*-**** 

1  TYPE  "INPUT  N. COMPLEX  X" 

ACCEPT  N.X 
PDT - 1 .0 
V=FLOAT(N) 

NuM=0 . 0 
DEN =0.0 
SCALE  =0 . 0 
AN= 1 . 0/X 
XINC=2.0/X 
COUNT =0 

C  IP  N  IS  CEPO,  USE  THE  INITIAL  VALLE  SIN(X)/S 

IF (N. EG . 01  GOTO  23 
1  7  4N  hN  rX  I  (L 

ituri=AN-Nun 

IFiNUM.EQ.0.0)  GOTO  IP 
PDT  =  PDT‘N'r- 
NUI  1  .  'J'  Nw:  ’ 

IF  CABS  i  REh_  i'PDT  j  .  .  lT.  1 .  D50)  GOTO  13 
P  DT  =P D  T  *  1 . D - 53 
SCHi_E=GCALE-SC 
10  COUNT “COUNT +1.3 

I F (COUNT . LE .VI  GOTO  19 
DEN=AN-DEN 

IF(DEN.EQ.C.O)  GOTO  19 

PDT=PDT/DEN 

DEN  = 1 .0/DEN 

IS  IF (COUNT.LT. 5.0)  GOTO  17 

C  THIS  CHECK  FOR  THE  COMPLEX  NUMERATOR  EQUAL  TO  THE  COMPLEX 

DENOMINATOR  MAY  BE  MACHINE  DEPENDENT 

IF (REAL(NUM) .NE. REAL (DEN) . OR. AIMAG(NUM) .NE. A I MAG (DEN) )  GOTO  17 
23  JO=CS IN (X) yX 

JN= JO/PDT 

URITE( 10. 100)  JN, SCALE. X.V.COUNT 

100  FORMAT ( 1H  , " JN* " . 2 (D20 . 10) . "  POUER  OF  10  SCALE* n . F6 . 0, /, 1H  ,"X*", 

X  2 f E  1 4 . 6 ) ,  ''  N  = 11 , F6 . 0 ,  ''  COUNT*". FB.0) 

GOTO  1 
END 
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